Timesaving Double-Grid Method for Real-Space Electronic-Structure Calculations 



> 



o 



I 

o 
o 



Tomoya Ono and Kikuji Hirose 

Department of Precision Science and Technology, Osaka University, Suita, Osaka 565-0871, Japan 

(February 1, 2008) 

We present a simple and efficient technique in ab initio electronic-structure calculation utilizing 
real-space double-grid with a high density of grid points in the vicinity of nuclei. This technique 
promises to greatly reduce the overhead for performing the integrals that involves non-local parts 
of pseudopotentials, with keeping a high degree of accuracy. Our procedure gives rise to no Pulay 
forces, unlike other real-space methods using adaptive coordinates. Moreover, we demonstrate the 
potential power of the method by calculating several properties of atoms and molecules. 

PACS numbers: 31.15.-p, 02.70.Bf, 36.40.-c, 71.15.Hx 



So far, a number of methods for ab initio electronic- 
structure calculations entirly in real space have been 
proposed. They have some advantages compared with 
the usual plane- wave approach. The first one is that 
boundary conditions are not constrained to be periodic, 
e.g., nonperiodic boundary conditions for molecules and 
a combination of periodic and nonperiodic boundary con- 
ditions for surfaces. Even more important is that a tech- 
nique utilizing real-space double-grid is available where 
many more grid points are put in the vicinity of nuclei, 
so that the integrals involving rapidly varying pseudopo- 
tentials inside the core regions of atoms can be calculated 
with a high degree of accuracy. This kind of integration 
over numerous sampling points, however, requires a large 
sum of computational effort. In this Letter, we present 
a quite simple and efficient double-grid technique which 
yields a drastic reduction of the computational cost with- 
out a loss of accuracy in the framework of the real-space 
finite-difference method. This technique can also be ap- 
plied to the plane-wave approach, if the integration over 
the core region is implemented in real space. 

The double-grid employed here consists of two sorts 
of uniform and equi-interval grid points, i.e., coarse- and 
dense-ones, depicted in Fig. |l]by the marks "x" and 
respectively. The dense-grid region enclosed by the circle 
is the core region of an atom that is taken to be large 
enough to contain the cutoff region of non-local parts 
of pseudopotentials. Throughout this paper we postu- 
late that wave-functions are defined and updated only 
on coarse-grid points, while pseudopotentials are strictly 
given on all dense-grid points in an analytically or nu- 
merically exact manner. 

Let us consider inner products between wave-functions 
ipix) and non-local parts of pseudopotentials v{x) (see 
Fig. H). For simplicity, the illustration is limited to 
the one-dimensional case hereafter. The values of wave- 
functions on coarse-grid points (o) are stored in com- 
puter memory, and the values on dense-grid points (•) are 
evaluated by interpolation from them. The well-known 
values of pseudopotentials both on coarse- and dense- 
grid points (o) are also shown schematically. Then, from 



Fig. ||(a) one can see that only the values on coarse-grid 
points are so inadequate that the inner products can not 
be accurately calculated; the errors are mainly due to 
the rapidly varying behavior of pseudopotentials. On 
the other hand. Fig. ||(b) indicates that the inner prod- 
ucts can be evaluated to great accuracy, if the number of 
dense-grid points is taken to be sufficiently large and also 
if the values of wave-functions on dense-grid points are 
properly interpolated from those on coarse-grid points 
@. 

There are several interpolation methods for wave- 
functions, among which the simplest is linear interpo- 
lation. In this case, the values of wave-functions ijjj = 
ip{xj) on dense-grid points interpolated from the 

values '^j = ip{X,j) on coarse- grid points Xj as 



h — [xj — Xj) 



{Xj+i - Xj) 



(1) 



where h is the grid spacing of coarse-grid points. The 
inner product is assumed to be accurately approximated 
by the discrete sum over dense-grid points, i.e.. 
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where d is the "diameter" of the core region, Rj is the 
atomic position, vj = v-'{xj) — v{xj — i?/), 2Ncore + 1 
{2nNcore+^) is the number of coarse- (dense-) grid points 
in the core region, hdens is the grid spacing of dense-grid 
points, and n = h/ hdens, i-e., n — 1 is the number of 
dense-grid points existing between adjacent coarse-grid 
ones. Now, substituting Eq.(P in the r.h.s. of Eq.(||), we 
have 
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As shown in Eq.(^), the r.h.s. of the inner product (H) 
has been replaced with the summation over coarse-grid 
points inside the core region, which produces only a mod- 
est overhead in the computational cost. It should be 
remarked that the weight factors Wj arising from the 
interpolation are independent of the wave- functions, but 
dependent only on the well-known values of pseudopoten- 
tials on dense-grid points. Thus, if once computing the 
factors Wj every molecular-dynamics time-step, we do 
not have to recalculate them throughout self-consistent 
iteration-steps. The extension of the above procedure 
to the cases of higher-order interpolations is straightfor- 
ward. 

Fourier interpolation is of great interest, since the term 
representing the position of the atom can be factorized 
in the expression of Wj as the structural phase factor 
exp{—ik^Rj). Indeed, the interpolated values of wave- 
functions ipj on dense-grid points xj, i.e., 

i>,= J2 (5a) 

with 

^k = -; J2 (5b) 

are substituted into the r.h.s. of Eq.(||) to give Eq.(||), 
where the weight factors in this case are 

Naar„ 

w'j= J2 ^fee^'^^'^e-^'^''^ (6a) 

with 

Wk = ^ Vje-'''^''^hdens- (6b) 

Here vj = v(xj). An advantage to be stressed is that 
the calculation of making the table of Wk in Eq. (|6b|) for 
each atom species, which requires a time-consuming com- 
putational effort because of the summation over dense- 
grid points, has only to be carried out once at the early 
stage of the entire job. At that time, the usage of the 
fast Fourier transforms (FFT) will considerably reduce 
an amount of computational operation. 

We now examine the performance of our method 
through calculations of several properties of atoms and 
molecules. Hereafter, we obey the nine-points finite 
difference formula (i.e., the formula with N=4) for 
the derivatives arising from the kinetic-energy opera- 
tor, imposing a nonperiodic boundary condition of van- 
ishing wave- functions. The dense-grid spacing is fixed 
at hdens = h/9. The electronic charge-density, the 
Hartree potential, and the exchange-correlation poten- 
tial are assumed to be described only on coarse-grid 



points. Exchange-correlation effects are treated using 
the local-spin-density approximation pof of the density- 
functional theory. The norm- conserving pseudopotential 
of Bachelet, Hamann, and Schliiter (BHS) is em- 
ployed in a separable non-local form p^ . 

The convergence of the total energy for the hydrogen 
atom as a function of the coarse-grid cutoff energy is 
presented in Fig. ^ According to Ref . M , we defined an 
equivalent energy cutoff -E™'""" [= (vr/ft)^ Ry] to be equal 
to that of the plane- wave method which uses a FFT grid 
with the same spacing as the present calculation. The 
position where the atom is located is at the center be- 
tween adjacent grid points. The hydrogen atom is one 
of the most difficult atoms to treat, owing to the rapid 
oscillation of its s-state non-local pseudopotential, and 
consequently the total energy as a function of the cutoff 
energy calculated without any interpolation sharply os- 
cillates. On the other hand, our prescription with inter- 
polations drastically improves the results; the obtained 
total energies converge rapidly and monotonously as the 
cutoff energy increases. 

Fig. ^ shows the total-energy variation as a function of 
the displacement of the oxygen atom relative to coarse- 
grid points along a coordinate axis. The coarse-grid spac- 
ing is taken toheh = 0.21 a.u. The first-row elements are 
difficult atoms to deal with, because their pseudopoten- 
tials are rapidly varying in the cutoff region. In treating 
these elements in the context with a real-space approach, 
it is necessary to consider the dependence of the energy 
on the position of the atom. As seen in Fig. ^, the en- 
ergy variation in our scheme with cubic interpolation is 
^^0.04% of the total energy, which is negligibly small. 

We next apply our method to the calculation of the 
equilibrium bond length of the CO2 molecule in order to 
examine the performance of our method for molecules. 
The total energies as a function of the C-0 distance are 
shown in Fig. y. We employ the cubic interpolation for- 
mula of the present method and take the coarse-grid spac- 
ing h to be 0.21 a.u. The results without any interpo- 
lation have many "humps" , and the distances between 
adjacent humps on the respective dotted curves are close 
to the grid spacing h, which confirms that the oscillation 
depends on the relative position of the atom with respect 
to coarse- grid points. To circumvent this problem, Gygi 
et al. [0 proposed a real-space grid in adaptive curvi- 
linear coordinates. However, owing to the use of these 
distorted coordinates, their scheme needs Pulay forces in 
molecular-dynamics simulations; it would be computa- 
tionally very demanding. On the contrary, our approach 
requires only the Hellmann-Feynman forces, which sig- 
nificantly reduces the computational cost concerning the 
calculation of forces. As shown in Fig. |^, the total ener- 
gies computed with our method do not make absurd de- 
pendence on the position of the atom. The bond length 
determined by the minimum of the total energy is in ex- 
cellent agreement with the experimental data 2.19 a.u. 
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These results make it clear that our method solves this 
problem and that it is very efficient and applicable . 

Finally, the calculated properties of the other 
molecules are given in Table |. The coarse-grid reso- 
lution h is 0.24 a.u. for N2 and 0.21 a.u. for the other 
molecules, and the cubic interpolation formula is used. 
Our results are in good agreement with those of experi- 
ments and other theories. 

In summary, we have presented a method for perform- 
ing the finite-difference electronic-structure calculations 
using real-space double-grid. Our method has the fol- 
lowing desirable features: (i)The inner products between 
wave- functions and pseudopotentials are evaluated with 
the same accuracy as in the case of dense-grid points, in 
spite of the calculation with respect to coarse-grid points. 
(ii)The computational effort is modest, thanks to the in- 
tegration over coarse-grid points. (iii)The double-grid 
acts as stabilizer in the calculation of the total energy, i.e., 
it suppresses the spurious oscillation of the energy for the 
displacement of the atom relative to coarse-grid points. 
(iv)Unlike other real-space methods using adaptive co- 
ordinates, our procedure gives rise to no Pulay forces, 
and so a substantial increase in the computational cost 
does not occur. From what has been mentioned above, 
it seems reasonable to conclude that our method is suit- 
able, because of its simplicity and efficiency, for large- 
scale molecular-dynamics simulations incorporating the 
norm-conserving pseudopotentials. Work is in progress 
to apply the method to large-scale molecular-dynamics 
simulations. 
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Solid State Physics at the University of Tokyo. Thanks 
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FIG. 1. Double-grid adopted in the text. The marks "x" 
and "•" correspond to coarse- and dense-grid points, respec- 
tively. The circle shows the core region of an atom which 
is taken to be large enough to contain the cutofl' region of 
non-local pseudopotentials. 

FIG. 2. Functions on coarse- and dense-grid points in the 
one- dimensional case. Xj (xj) represents a coarse- (dense-) 
grid point with j — nj -\- s(0 < s < n), and so Xj = Xnj. (a) 
The inner product between wave-function ij^ix) and non-local 
part of pseudopotential v{x) evaluated over coarse-grid points, 
(b) The inner product evaluated over dense-grid points. 

FIG. 3. Convergence of the total energy for the hydro- 
gen atom as a function of the coarse-grid cutoff energy 
^coars 'YYie atom is located at the center between adjacent 
coarse-grid points. 
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FIG. 4. Variation of the total energy as a function of the 
displacement of the oxygen atom relative to coarse-grid points 
along a coordinate axis. The coarse-grid spacing h is 0.21 a.u. 
The result with linear interpolation is coincident with that of 
the cubic case within the relative error of 0.5%. 

FIG. 5. Total energy of the CO2 molecule in the sym- 
metric linear conformation and location of the carbon atom. 
(a)The carbon atom is located at the center between adjacent 
coarse-grid planes in each direction. The marks "x" repre- 
sent coarse-grid points. (b)The carbon atom is located on 
the coarse-grid plane which is at right angles to the bonding 
axis and at the center between adjacent coarse-grid planes in 
the other directions. (c)Total energy of the CO2 molecule as 
a function of the C-0 distance. The coarse-grid spacing h 
of 0.21 a.u. is the same size as the corresponding Euclidean 
grid spacing in the adaptive coordinate approach The 
calculated points are fit to spline functions as a guide to the 
eye. 



TABLE I. Properties of diatomic molecules. The experi- 
mental data are from The theoretical results are from 
our present calculation and from other methods using similar 
forms for the local-spin-density approximation. 





N2 


O2 


CO 


Bond length (a.u.) 








Experiment 


2.07 


2.28 


2.13 


This work 


2.07 


2.28 


2.13 


Other theory 


2.07*" 


2.27*^ 


2.13^^ 


Vibrational mode (cm~^) 








Experiment 


2358 


1580 


2169 


This work 


2390 


1500 


2200 


Other theory 


2380'' 


1620'" 


215r 


Cohesive energy (eV) 








Experiment 


9.76 


5.12 


11.09 


This work 


11.20 


7.57 


12.80 


Other theory 


11.6*^ 


7.6'^ 


12.8'' 



^From Ref. [Rk 
''From Ref. 
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